Quasi-Spherical Approximation for Rotating Black Holes 
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We numerically implement a quasi-spherical approximation scheme for computing gravitational 
waveforms for coalescing black holes, testing it against angular momentum by applying it to Kerr 
black holes. As error measures, we take the conformal strain and specific energy due to spurious 
gravitational radiation. The strain is found to be monotonic rather than wavelike. The specific 
energy is found to be at least an order of magnitude smaller than the 1% level expected from typical 
black-hole collisions, for angular momentum up to at least 70% of the maximum, for an initial 
surface as close as r = 3m. 
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I. INTRODUCTION 

A quasi-spherical approximation scheme in a 2+2 de- 
composition of the space-time has recently been intro- 
duced 0. This proposal is with the aim of providing 
a computationally inexpensive estimate of the gravita- 
tional waveforms produced by a black-hole or neutron- 
star collision, given a full numerical simulation up to (or 
close to) coalescence, or an analytical model thereof. 

The scheme truncates the Einstein equations by remov- 
ing second-order terms which would vanish in a spheri- 
cally symmetric space-time, cf. Bishop et al. Thus 
when the linearized fields vanish, spherical symmetry is 
recovered in full. Unlike previous work on null-temporal 
formulations ||,[||, a dual- null formulation is adopted 
here, i.e. a decomposition of the space-time by two inter- 
secting foliations of null hypersurfaces. Technical advan- 
tages of the scheme include that only ordinary differential 
equations need be solved, and that the dual-null formu- 
lation is adapted to radiation extraction. Advantages 
of applicability include that no prescribed background is 
required and that arbitrarily rapid dynamical processes 
(close to spherical symmetry) are allowed. The pressing 
question concerns how well the scheme handles deviations 
from spherical symmetry. The principal such deviation 
in the context of coalescing black holes is expected to 
be due to angular momentum. The primary test case is 
therefore Kerr black holes, the unique stationary vacuum 
black holes. 

This article reports a numerical implementation of the 
quasi-spherical approximation and its application to Kerr 
black holes, taking Boyer-Lindquist quasi-spheres. Ex- 
cept in the non-rotating (Schwarzschild) case, the ap- 
proximation produces spurious gravitational radiation, 
not present in the exact solution. We consider two 
measures of the error introduced by the approximation. 
Firstly, the practical measure is the waveform of the spu- 
rious strain, as compared to signals expected to be mea- 
sured by interferometers. Secondly, we measure the spe- 



cific energy, i.e. the ratio of the radiated energy E to the 
original mass m; this is a conservative measure, as it in- 
volves summing the errors over all angles of the sphere, 
whereas observation is restricted to a particular angle. 

In principle these quantities depend on only the spin 
parameter a/m, the relative initial radius ro/m and, in 
the case of the strain, the angle. For the approximation 
to be useful, the error should be significantly less than the 
values expected for a realistic black-hole collision. Typi- 
cal values for the specific energy obtained from numerical 
simulations Q] or from the close-limit approximation (|] 
have increased from early estimates to around 1% if the 
initial relative momentum Q or angular momentum JtJ 
is appreciable. The theoretical limit on how badly an ap- 
proximation might perform is much higher: 29% of the 
mass of a maximally rotating Kerr black hole may be 
extracted by the Penrose process [|| . 

The article is organized as follows. Section || describes 
the dual-null formalism, the quasi-spherical approxima- 
tion and the observables, strain and energy. Section 



[II 



describes our model and numerical integration pro- 
cedures. The numerical results are shown in Section IV 
and we summarize the article in Section [v| 



II. FORMULATION AND APPROXIMATION 

A. Dual-null formulation 

The quasi-spherical approximation jl| is based on a 
dual- null formulation H of Einstein gravity [jl0|, sum- 
marized as follows. One takes two intersecting families 
of null hypersurfaces labelled by x^. Then the normal 
1-forms n = —dx^ satisfy 
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(2.1) 



where g is the space-time metric. The relative normal- 
ization of the null normals may be encoded in a function 
/ defined by 
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e / = -p- 1 (n+,n"). (2.2) 

Then the induced metric on the transverse surfaces, the 
spatial surfaces of intersection, is found to be 

h = g + 2er f n + <g> n~ , (2.3) 

where (g> denotes the symmetric tensor product. The co- 
variant derivative of h is denoted by D. The dynamics is 
described by Lie transport along two commuting evolu- 
tion vectors u±: 

[u+,u_]=0. (2.4) 

Specifically, the evolution derivatives, to be discretized 
in a numerical code, are 

A±=1_L U± , (2.5) 

where _L indicates projection by h and L denotes the Lie 
derivative. There are two shift vectors 

s± = ±u±. (2.6) 

In a coordinate basis (u+, tt— ; e,-) such that u± = d/dx , 
where = d/dx 1 is a basis for the transverse surfaces, 
the metric takes the form 

g = hij(dx l + s\dx + + s l _dx~) ® 

(dar* + s J + dx + + s 3 _dx~) - 2e~ f dx + ® dx~. (2.7) 

Then (/i, /, s±) are configuration fields and the indepen- 
dent momentum fields are found to be linear combina- 
tions of 

6± = *L ± *l, (2.8) 
a± = ±Z±h-6±h, (2.9) 
v± = L±f, (2.10) 
w = ±efh([l-,l+]), (2.11) 

where * is the Hodge operator of h and L± is shorthand 
for the Lie derivative along the null normal vectors 

l± = u± - s± = e- f g- 1 (n T ). (2.12) 

Then the functions 9± are the expansions, the traceless 
bilinear forms <j± are the shears, the 1-form to is the 
twist, measuring the lack of integrability of the normal 
space, and the functions v± are the inaffinities, measur- 
ing the failure of the null normals to be affine. The fields 
(0± ,cr±,u±,LL)) encode the extrinsic curvature of the dual- 
null foliation. These extrinsic fields are unique up to du- 
ality ± i — s- =p and diffeomorphisms which relabel the null 
hypersurfaces, i.e. dx^ i— ► e x± dx for functions X±(x ± ). 

It is also useful to decompose h into a conformal factor 
n and a conformal metric k by 

h = Q- 2 k, (2.13) 

such that 



A±*1 = 0, (2.14) 

where * is the Hodge operator of k, satisfying *1 = *Q~ 2 . 
Taking quasi-spherical coordinates x 1 — (0, 0) such that 
il = sin 9d0 A d(f>, SI -1 is the quasi-spherical radius. In 
an asymptotically flat space-time, it becomes convenient 
to use the conformally rescaled expansions and shears 

d ± = Q- 1 9±, (2.15) 

s± = n<j±, (2.16) 

since they are finite and generally non-zero at null infinity 

B. Quasi-spherical approximation 

Of the dynamical fields and operators introduced 
above, (s±, a±, to, D) vanish in spherical symmetry, while 
(h, /, 0±, v±, A±) generally do not. The quasi-spherical 
approximation consists of linearizing in (s±, a±, w, D), 
i.e. setting to zero any second-order terms in these quan- 
tities. This yields a greatly simplified truncation ]IJ] of 
the full field equations, the first-order dual-null form of 
the vacuum Einstein system |Tc| ]. In particular, the trun- 
cated equations decouple into a three-level hierarchy, the 
last level being irrelevant to determining the gravita- 
tional waveforms. The remaining equations are the quasi- 



spherical equations 

A±ft = -ifi 2 i?±, (2.17) 

A±/ = i/±, (2.18) 

A±z?± = -v±<&±, (2.19) 

A±tf T = -ft(iiM_ + e _/ ), (2.20) 

A±v T = -ft 2 (±tf+tf_ + e~t), (2.21) 

and the linearized equations 

A±fc = n<f±, (2.22) 

A±? T = fi(?+ ■ k- 1 • ? _ - |^?±). (2.23) 



These are all ordinary differential equations; no trans- 
verse D derivatives occur. Thus we have an effectively 
two-dimensional system to be integrated independently 
at each angle of the sphere. 

The initial-data formulation is based on a spatial sur- 
face S orthogonal to l± and the null hypersurfaces S± 
generated from S by l±, assumed future-pointing. The 
initial data for the above equations are (Q, /, k, &±) on 
S and ( i z±,v±) on £±. We will take 1+ and /_ to be 
outgoing and ingoing respectively. 



C. Strain 

The variables are directly related to physically mea- 
surable quantities. In particular, 
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is the transverse strain tensor measured along a worldlinc 
7 normal to the transverse surfaces, where k = -LL\k in 
terms of a vector A = d/dr tangent to 7. For a detec- 
tor at large distance, one may apply the linearized ap- 
proximation, where 2e reduces to the transverse traceless 
metric perturbation of a linearized plane gravitational 
wave. In a weak gravitational field, one may use New- 
tonian physics, where e reduces to the Newtonian strain 
tensor. Thus the displacements to be measured by an 
interferometer are 



T = e ( e ' e ) 



(2.25) 



where the unit vector e is the direction of displacement. 
Writing A = a + l + + a~Z_ yields 



f 



f2(a + Cf + a <^_)dr. 



(2.26) 



Since the strain vanishes at future null infinity S + , it is 
convenient to use the conformal strain tensor 



£ — — I (I.r 



(2.27) 



where the integal is at constant x + . We will denote its 
plus and cross components by e+ = egg and e x = eg^. In 
order to compare with observational results, one converts 
back to the strain, 



(2.28) 



where R is the distance between the source and the de- 
tector. 



intersection with £+, which in the Kerr case will be just 
the ADM mass m. We propose using the specific energy 
E/m at 3 + as a measure of the strength of a black-hole 
collision, following various references . This is the 
fraction of the original mass-energy which has radiated 
away. 



III. MODEL AND NUMERICAL PROCEDURES 



A. Model: Kerr black hole 



As our model, we take a Kerr black-hole geometry, 



ds 2 = ~ — [dt — asin 2 aj - |J 

2-i 



A 



dr 2 + Y,d6\ 



where 



A = r — 2mr + a 
£ = r 2 + a 2 cos 2 ' 



-[(r 2 + a 2 )d<P - adt} 2 
(3.1) 



(3.2) 
(3.3) 



and m is the mass and am the angular momentum. The 
horizon radius is denoted by rjj = m + V m 2 — a 2 . We 
take the quasi-spherical approximation adapted to these 
(Boyer-Lindquist) coordinates, i.e. the initial surface S is 
of constant r — r and constant t, as depicted in Fig^. 
The inaccuracy of the quasi-spherical approximation as 
measured by E/m then depends in principle only on a/m 
and ro/m. We expect E/m to be monotonically increas- 
ing in a/m, from zero at a = to a maximum at a = m, 
since angular momentum is the cause of the asphericity. 
Similarly, E/m is expected to be monotonically decreas- 
ing in ro/m, to zero at infinity, since the approximation 
should be better at large distances. 



D. Energy 



B. Initial data 



We define the energy flux <f> of the gravitational waves, 
or more conveniently, the conformal energy flux tp — 
fl~ 2 (j), as the 1-form tp = ip + dx + + ip-dx~ , where 



e f d^k ah k cd 
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S±acS±bd 



64tt 



(2.29) 



These expressions have the same form as those for the 
conformal Bondi flux at 3^ fjj], but we propose using 
them locally. Then 0_ is the outgoing flux and <fi+ is the 
ingoing flux. The corresponding energy E of the gravi- 
tational waves is then given by 



A±E 



*<P± 



(2.30) 



with the initial condition E\g = 0. Thus the Bondi en- 
ergy at S + is E + E , where E is the Bondi energy at the 



An explicit dual-null form of the Kerr metric is not 
known except on the symmetry axis ]TT[ ] or in the 
Schwarzschild the Kruskal form. Although there 

is an effort for this direction Jl2] ), wc did not find ex- 
plicit double-null coordinates which arc well behaved at 
the outer horizons and infinity, despite trying changes 
in angular coordinate and Kruskal-type rescalings p3[ . 
However, we can construct the initial data analytically 
as functions of (r, 6) , then convert to the required func- 
tions of {x ± 1 9) 1 as follows. We remark that our initial 
surfaces are £±, and the following method applies out- 
side the horizons. 

The null normal vectors are initially given by 



2(A- 



a 2 sin 2 t 



1/2 f A \ 1/2 



3 



r = r H/ 



r = Xr 




nr H 



FIG. 1. The region of numerical integration is shown as 
the shaded region in the picture. Initial data is prescribed on 
a spatial surface S of constant Boyer-Lindquist r — ro and 
t, and the null hypersurfaces E± generated from it. On £_ 
(£+), the x~ (x + ) coordinate is set so as to cover the region 
Xth < r < ro (ro < r < urn), where 1 < A ft* 1 and n>l 
are constants to be set by hand. 



where the normalization is such that 



/I, 



0. 



(3.5) 



The apparent degeneracy at A = a 2 sin 2 9 is just the 
boundary of the ergoregion where dt becomes spatial; 
the dual-null coordinates extend through. We also fix 



0. 



(3.6) 



which means that x |s ± are affine parameters. This im- 
plies /|s ± = 0, fixing Zzp|s ± and therefore locally deter- 
mining the dual-null foliation. 

The quasi-spherical conformal factor is 



n = ((r 2 + a 2 ) 2 - Aa 2 sin 2 9)^ 



(3.7) 



which is real and positive. Then we obtain the conformal 
metric 



k = VL 2 Hd9 2 + 



sin 



(3.8) 



and the conformally rescaled expansions and shears 
"A 




n 3 u, 

Q 5 Va 2 sm 2 9 f d9 2 



sin 2 9 
f^S 2 ' 



(3.9) 
(3.10) 



where 

U = -2fl- 5 d r n = 2r(r 2 + a 2 ) - (r - m)a 2 sin 2 9, (3.11) 
rr 6 <9 r (f7 2 £) 



V 



a 2 sin 2 i 



+ 3mr 2 + a 2 (r-m) cos 2 9. (3.12) 



To complete the initial data construction, we need to 
know the initial data on S± as functions of [x , 9). So 
we need to know r|s ± as functions of (x ,6). This is 
determined by the equations 



A±r| s± =± 



_A 

2E 



1/2 



giving 



= ± 



(3.13) 



(3.14) 



where the integral is along a curve of constant (6,<f>). 
Note that the A factor means that we take S outside 
the horizons r = rjj, which anywa y is t he region of in- 
terest. We numerically integrated (3.14) using a fourth 
order Runge-Kutta method (Fehlberg method), then in- 
verted. This was checked against the analytic solution in 
the equatorial plane 9 — it/ 2: 



a± = ±V2 (VA + roln(r 



VA)Wc. (3.15) 



As to is an overall scale, we fixed it to unity. 



C. Evolution procedures 

Here we describe our numerical procedures. We have a 
set of ordinary differential equations in two variables. As 
pointed out by Gundlach and Pullin |l4|| , free evolution 
schemes in such a system may lead to unstable evolution. 
This fact was also seen in our experience, and we devel- 
oped a kind of predictor-corrector scheme similar to that 
of Hamade and Stewart (lal . 




FIG. 2. The dual-null integration scheme. In order to ob- 
tain the data at grid u^ +1 , we need both u£ and u^j. 

The actual steps we took are the followings. The set of 
variables is u = (O, /, t?±, v±, k a b, <T±ab)- Let us schemat- 
ically express a set u at a point x~ = A; on a slice 
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= n) as Uu- The data u^ +1 is determined from 
both and ujj as in Fig. ^. Suppose we have already 
all the data at an d uJJ- 

(1) Firstly, we evolve along the x + -direction, say from 
M fc t° We have a set of equations for 



A+ft = 




(3.16) 


A+/ = 




(3.17) 


A+^+ = 




(3.18) 


A + i?_ = 




(3.19) 


A+i/_ = 


-fi 2 (§0+0_ +e-'), 


(3.20) 


A+fc ab = 




(3.21) 


A + <r_ ah = 




(3.22) 



results shown in this article are obtained by setting the 
tolerance in the above step (3) to 10 -5 . The code was 
tested for the Schwarzschild case, for which the analytic 
expression in dual-null coordinates is known; the calcu- 
lated expansion &± differed from the exact expression to 
within 10~ 6 . 

In the next section, we present our evolutions of a Kerr 
black-hole space-time under this quasi-spherical approx- 
imation. We chose the initial null slice £_ so as to cover 
the region 1.25r# < r < r . We stopped the evolution 
at x + — 30, which corresponds to r being 25 ~ 30m, de- 
pending on the values of tq and a. We took 51 grid points 
in the x~ direction and 11 grid points in 9 = [0, 7r/2], and 
evolved with grid separation Ax + — 0.5Ax~. 

IV. NUMERICAL RESULTS 



The step is integrated using the Fehlberg method. 
Note that we do not have equations for evolving v+ 
and Cf j therefore we have to interpolate them using 
{v+^+ab)k and { v +^+ab)k +1 ■ Tne latter was lin- 
early extrapolated for the first iteration, but will be 
updated after an integration along the a; _ -direction 
(next step) has been done. 

(2) Secondly, we evolve along the ir~-direction, from 
We have a set of equations for 



,n+l 

1 k-i 



(3.23) 



to Ufe 

A_i/ + = -ft 2 (±i?+i?_ +e-f) 
A^ +ab = n(s +ac k cd ^ db - §tf+?-ab), (3.24) 

for completing the set u, but we also evolve #±, f2 
and / by 



A_ft 
A-/ 

A_tf_ 
A_i?+ 



(3.25) 
(3.26) 
(3.27) 
(3.28) 



Here again we have to interpolate and v_ in in- 
tegrating the above, and we use a cubic spline inter- 
polation using S- a b)k +1 (1 < h < k), where 



the data (z/_, <r_ a6 )£ H 
step (1). 



was given in the previous 



monitoring the differences of (k ab , i?±, f)% from 
the above step (1) and (2). If they are all within 
a tolerance, then we finish this evolution step by 
updating (y + , s+ a b, #±> /) as a value at 
not, we repeat back to the step (1). 



,n+l 



If 



We construct a numerical grid in x^ space with con- 
stant spacing in each direction. The iteration procedures 
are completed a couple of times at each grid point. The 



(3) We check the consistencies of the evolution, by (^ 



Recall that our principal measures of the gravitational 
radiation are the co nform al strain e (2.27) and the spe- 
cific energy E/m C^TaoT) . Since we are using confor- 
mal variables, we expect that we can evolve towards the 
asymptotically flat region without long-term evolution in 
the x + direction. In Figj3[ we plotted the specific energy 
E/m at the boundaries of the i ntegr ation region. We 
integrated the A + E equation of ( [2.30 ) along the hyper- 
surface S + (x~ — 0), setting E = on S (x + = x~ =0), 
then integrated E using the A_E equation at each con- 
stant x + . We plotted E as a function of x~ at a constant 
x + surface in Fig.||(a). We see that E is converging to a 
line (the solid line in the figure), and not diverging even 
close to the black hole (at larger x' -). Fig.|(b) plots E 
at the final value of x~ as a function of x + . We see from 
the figure that the energy measured for increasing x + 
converges at some value, as expected. 

For the same set of parameters, we also plot the evolu- 
tion behaviour of the conformal strain in Fig.^ and Fig||. 
The cross component, e x = £g0, is zero in this model, so 
only the plus component e + = eg g is n eeded. The con- 
formal strain is calculated from ( |2.27 ) as a function of 
x~ at constant x + by setting e + = at x~ — 0. We 
again observe that e+ converges to particular lines (the 
solid lines) as x + increases, again reflecting the confor- 
mal variables. The line of x + — 30.0 in Fig.[|, therefore, 
is close to the waveform for observers infinitely far from 
the source. 

Note that the horizontal axis in Fig.[|is x~ coordinate, 
and this is related to Ax~ w y/2At at large distance, cf. 

15). Then the unit length Aa; - for observers at large 
distance is about 5 (m/M Q )/i sec, translated from our 
units c = G = 1. Our plot, therefore, covers a quite 
short time period compared with the typical millisecond 
timescale of gravitational waves from a Kerr black hole.^] 



'According to the quasi-normal mode analysis of the Kerr 
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FIG. 3. Specific energy E/m for a/m = 0.1. (a) E/m 
is plotted as a function of the ingoing null coordinate x~ for 
each constant outgoing null coordinate x + — 0.0, 3.0, • • • , 30.0. 
We set ro = 4.0 for this plot, (b) The integrated E/m over 
the ingoing null coordinate x~ is shown as a function of the 
outgoing null coordinate x + . We show both ro = 3.0 and 4.0 
cases. We see that the specific energy converges to a particular 
positive value in the x + direction, as expected. 



To obtain longer timescales we would have to integrate 
closer to the horizon, which causes numerical difficulties 
due to the infinite redshift. However, for a dynamically 
evolving black hole, the event horizon has finite redshift 
and so could lie in the numerical integration region, al- 
lowing evolution to late times. 

The magnitude of the conformal strain in Fig.^j is 
rescaled to the observable strain by (2.28). We can 



black hole [[I6| , the dominant frequencies (fundamental mode 
corresponding to / = 2) of quasi-normal mode for a 10 Mq 
black hole is between 1.2 kHz (for a = 0) and 1.8 kHz (for 
close to maximally rotating). 



FIG. 4. Conformal strain e+ for a/m = 0.1 and ro = 4.0. 
The plot is for the equatorial plane 9 = 7r/2, showing the 
convergence of these lines in the x + direction. Lines are of 
x + = 3.0, 6.0, 9.0, • ■ ■ , and 30.0. We remark that these lines 
are not wave- like. 



compare with an example of expected strain e ~ 10 20 
for R = 100 Mpc by converting our units: e = 
3 x 10 _27 e/(i?/100 Mpc). This is small enough to vali- 
date the quasi-spherical approximation. The converged 
conformal strain (x + = 30 lines in Fig.Q) is increasing 
as the ingoing coordinate x~ approaches the black hole 
horizon. However, if we extrapolate this magnitude to 
the horizon (which will be reached around x~ ~ 6.0 for 
this choice of parameter), it is still many orders of mag- 
nitude less than expected values. 

Moreover, the strain does not behave like a wave in 
this case. This is good news for future applications of 
the quasi-spherical approximation, because the produced 
spurious waveform is quite different from a normal grav- 
itational wave. We also show 6 and a dependencies of e+ 
in Fig.| 

Our final, most conservative check of the quasi- 
spherical approximation is to compare the specific energy 
E/m with the expected specific energy of gravitational 
waves from an inspiralling black-hole binary. The Kerr 
black-hole space-time seems to be a good example for 
comparing with the result of the close-limit approach [Jt) . 
In Fig.[| we plotted the specific energy E/m due to spuri- 
ous radiation, as a function of a/m and ro /m. We applied 
the same grid points and other parameters in numerics 
with previous figures, and evaluated E/m at x + = 30. 
For higher a and larger ro cases, we could not fill plots 
in FigJq. This is because we kept the resolutions and the 
same tolerance for the consistency convergence criteria 
for all cases, and these criteria failed for higher a and 
larger ro. If we increase the resolutions and/or adjust 
the convergence criteria, then we can fill in these missing 
points also. 

Consequently, we observe that the specific energy E/m 
increases with a/m and decreases with ro/m, as ex- 



G 



FIG. 6. Logarithmic plot of specific energy E/m due to 
spurious radiation, as a function of a/m and ro/m. En- 
ergy is measured at x + — 30, and the plotted range is 
ro/m £ [3.0,4.5] and a/m G [0.1,0.7]. 




FIG. 5. Conformal strain e+ for the same parameters as 
Fig.|. (a) is e+ at x + = 30.0 for different 8. We see that the 
maximal strain occurs in the equatorial plane, 9 = 7i"/2, as 
expected, (b) shows the dependence of e+ on a/m. The lines 
are for the data at x + — 30.0 for 6 — tt/2. Both solid lines 
are equivalent with the solid line in FigMl (a). 



pected. If we compare the amplitude of E/m with Fig.l 
of Khanna et al. fft, then we find that our values are 
at least an order of magnitude smaller than the results 
of the close-limit approximation, up to the range where 
different versions of the latter diverge. The fact that 
the spurious radiation produced in the quasi-spherical 
approximation is quite small indicates the robustness of 
this approximation to the general situations. 



V. CONCLUDING REMARKS 

We tested the quasi-spherical approximation by apply- 
ing it to Kerr black holes. We numerically calculated the 
strain and energy flux of the spurious gravitational ra- 
diation produced from this approximation, and showed 
that (a) it converges quickly due to our conformal vari- 



ables, (b) it does not behave like wave-like oscillations, 
and (c) the total radiated energy is at least an order of 
magnitude less than the gravitational radiation emission 
estimated from coalescing binary black holes, according 
to the close-limit approximation 0. We remark that 
the close-limit approximation is the only current result 
which predicts the total amount of radiation from inspi- 
ralling binary black holes. Numerical results for head-on 
collisions with appreciable relative momentum also give 
similar estimates 

These results suggest that the spurious radiation does 
not fatally affect the gravitational waveform estimation. 
It might not affect the waveform estimation at all, and we 
might extract its effect to the total energy by subtract- 
ing the amount we showed in Fig.[| These facts directly 
encourage the robustness of the quasi-spherical approx- 
imation. Therefore we are interested in applying this 
scheme to more general situations, and/or implementing 
it as an output routine for full numerical simulation codes 
of binary black holes or compact stars, such as those us- 
ing the standard 3+1 decomposition of spacetime. These 
efforts will be reported elsewhere. 
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